Slotted Diabolo Stripline

# standard python imports
import os
import numpy as np
import gdstk
import time

# # plot libraries
import matplotlib.pyplot as plt
import plotly.graph_objects as go
import plotly.io as pio
import plotly.express as px
# pio.templates.default = "simple_white"
pio.templates.default = "plotly_dark"
pio.renderers.default = "notebook"

# tidy3d import
import tidy3d as td
from tidy3d import web

# users packages
from metasruface import make_wings, make_wings_array, make_diabolo, make_diabolo_array
from utility import plot_gdscell
from material.THzMaterial import diamond_thz_nd, silica_thz_nd, pdms10_thz_nd

FLOAT_MIN = 4.9406564584124654E-324

Draw the Metasurface in GDS

ncol_ms = 11 # array size
nrow_ms = 11
d_ms = 25.0 # bridge width
l_ms = 135.0 # unit size
p_ms = 270.0 # periodicity 
s_ms = -25.0 # separation between two triangles
theta_ms = 50.0 # flare angle [deg]
l_magbox = d_ms/np.tan(theta_ms/2.0/180*np.pi)+s_ms
ld_thz_ms = ld_thz_ms= {"layer": 11, "datatype": 1}

cell = make_wings(d_ms, l_ms, p_ms, s_ms, theta_ms, **ld_thz_ms)
cellcell = make_wings_array(ncol_ms, nrow_ms, d_ms, l_ms, p_ms, s_ms, theta_ms, **ld_thz_ms)

um_gds = 1E-6
gdslib = gdstk.Library()
gdslib.unit = um_gds
gdslib.add(cellcell, *cellcell.dependencies(True))
filename = f"temp_{20240802}.gds"
gdslib.write_gds(f"gds_output/{filename}")

# Load a GDSII library from the file we just created
gds_path = f"gds_output/{filename}"
lib_loaded = gdstk.read_gds(gds_path)

# Create a cell dictionary with all the cells in the file
all_cells = {c.name: c for c in lib_loaded.cells}
print("Cell names: " + ", ".join(all_cells.keys()))
cell_loaded = all_cells["Wings"]
print(cell_loaded)
plot_gdscell(cell_loaded)

Cell names: Diabolo Array, Wings
Cell 'Wings' with 2 polygons, 0 flexpaths, 0 robustpaths, 0 references, and 0 labels

Simulation Setting and Visualization

# define freqeunecy range
GHz = 1e9  # 1 GHz = 1e9 Hz
freq0 = 400.0 * GHz  # central frequency
freqs = np.linspace(freq0-200.0*GHz, freq0+400.0* GHz, 500)  # frequency range of interest
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))  # width of the source spectrum
lda0 = td.C_0 / freq0  # central wavelength


# define structure geometric parameters
t_sub = 180.0 #um
t_diabolo = 0.15
t_magbox = 8.0
t_magbox_open = 3.5
t_rfstrip = 1.0
w_rfstrip = 100.0
t_uf = 60.0
t_ufch = 30.0

# define material properties, assume non-dispersive dielectric
sigma_copper = 42  # copper conductivity in S/um, ref DOI:10.1109/TTHZ.2015.2468074
t_copper = t_diabolo
copper3D = td.Medium(conductivity=sigma_copper)
copper2D_thz = td.Medium2D.from_medium(copper3D, thickness=t_diabolo)  # define copper as a Medium2D
copper2D_rf = td.Medium2D.from_medium(copper3D, thickness=t_rfstrip)  # define copper as a Medium2D

sub_medium = silica_thz_nd  # define substrate medium
msbox_medium = diamond_thz_nd  # define medium in MS box
uf_medium = pdms10_thz_nd  # define microfluidic medium

reference_plane = "bottom" # Reference plane where the cross section of the device is defined
z_fp_compen = 0.000001
msatom_geo = td.Geometry.from_gds(
    cell_loaded,
    gds_layer=ld_thz_ms["layer"],
    gds_dtype=ld_thz_ms["datatype"],
    axis=2,
    # slab_bounds=(0.0, t_diabolo),
    slab_bounds=(z_fp_compen, z_fp_compen), # assume thin structure, the effect of finite thickness is handled in Medium2D
    reference_plane=reference_plane,
)

rfstrip_geo = td.Box(
    center=(0, 0, t_ufch+t_rfstrip/2.0), 
    size = (w_rfstrip, td.inf, 0)
) #assum 2D metal to reduce simulation time

# msatom = td.Structure(geometry=msatom_geo, medium=copper3D)
msatom_wings = td.Structure(geometry=msatom_geo, medium=copper2D_thz, name="MS Wings")


msatom_flat = td.Structure(
    geometry=td.Box(center=(0, 0, t_magbox), size=(l_magbox, d_ms, 0)),
    medium=copper2D_thz,
    name="MS Flat"
)

msatom_wingside_1 = td.Structure(
    geometry=td.Box(center=(l_magbox/2.0, 0, t_magbox/2.0), size=(0, d_ms, t_magbox)),
    medium=copper2D_thz,
    name="MS Wing SideWall +x"
)

msatom_wingside_2 = td.Structure(
    geometry=td.Box(center=(-l_magbox/2.0, 0, t_magbox/2.0), size=(0, d_ms, t_magbox)),
    medium=copper2D_thz,
    name="MS Wing SideWall -x"
)

msatom_band_1 = td.Structure(
    geometry=td.Box(center=(0, d_ms/2.0, t_magbox/2.0+t_magbox_open/2.0), size=(l_magbox, 0, t_magbox-t_magbox_open)),
    medium=copper2D_thz,
    name="MS Band +y"
)

msatom_band_2 = td.Structure(
    geometry=td.Box(center=(0, -d_ms/2.0, t_magbox/2.0+t_magbox_open/2.0), size=(l_magbox, 0, t_magbox-t_magbox_open)),
    medium=copper2D_thz,
    name="MS Band -y"
)

msatom = [msatom_wings, msatom_flat, msatom_wingside_1, msatom_wingside_2, msatom_band_1 ,msatom_band_2]

rfstrip = td.Structure(geometry=rfstrip_geo, medium=copper2D_rf, name="RF Stripline")

substrate = td.Structure(
    geometry=td.Box(center=(0, 0, -t_sub / 2), size=(td.inf, td.inf, t_sub)),
    medium=sub_medium,
    name="Base Substrate"
)

magbox = td.Structure(
    geometry=td.Box(center=(0, 0, t_magbox / 2), size=(l_magbox, d_ms, t_magbox)),
    medium=msbox_medium,
    name="Mag Box"
)

ufluidics = td.Structure(
    geometry=td.Box(center=(0, 0, t_uf/2.0), size=(td.inf, td.inf, t_uf)),
    medium=uf_medium,
    name="Microfluidics"
)

offset = lda0 / 2  # extra spacing added in the positive and negative z directions
# define a plane wave source
plane_wave = td.PlaneWave(
    center=(0, 0, t_uf+0.1*offset),
    size=(td.inf, td.inf, 0),
    source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
    direction="-",
)

# define a flux monitor to measure reflection
S11_monitor = td.FluxMonitor(
    center=(0, 0, t_uf+offset),
    size=(td.inf, td.inf, 0),
    freqs=freqs,
    name="S11",
)

# define a flux monitor to measure reflection
S21_monitor = td.FluxMonitor(
    center=(0, 0, -t_sub - offset),
    size=(td.inf, td.inf, 0),
    freqs=freqs,
    name="S21",
    normal_dir="-",
)

# define a field monitor to visualize field distribution
field_monitor_xc0 = td.FieldMonitor(
    center=(0, 0, 0.0), size=(0, td.inf, td.inf), freqs=[freq0], name="field_xc0"
)

field_monitor_yc0 = td.FieldMonitor(
    center=(0, 0, 0.0), size=(td.inf, 0, td.inf), freqs=[freq0], name="field_yc0"
)

# define a field monitor to visualize field distribution
field_monitor_zc0 = td.FieldMonitor(
    center=(0, 0, 0.0), size=(td.inf, td.inf, 0), freqs=[freq0], name="field_zc0"
)

# define a field monitor to visualize field distribution
field_monitor_zct = td.FieldMonitor(
    center=(0, 0, t_magbox), size=(td.inf, td.inf, 0), freqs=[freq0], name="field_zct"
)

# define a field monitor to visualize field distribution
field_monitor_ms_boxcenter = td.FieldMonitor(
    center=(0, 0, t_magbox), size=(td.inf, td.inf, 0), freqs=[freq0], name="field_boxcenter"
)

# define a field monitor to visualize field distribution
t_roi = t_magbox*1.5
field_monitor_roi = td.FieldMonitor(
    center=(0, 0, t_magbox/2.0-(t_roi-t_magbox)/2.0*0.9), size=(d_ms, d_ms, t_roi), freqs=freqs, name="field_roi"
)

# define a BoundarySpec
boundary_spec = td.BoundarySpec(
    x=td.Boundary.periodic(),
    y=td.Boundary.periodic(),
    z=td.Boundary(minus=td.PML(), plus=td.PML()),
)

# define a mesh override structure
refine_ms = td.MeshOverrideStructure(
    geometry=td.Box(center=(0.0, 0.0, 0.0), size=(l_ms*1.2, (l_ms-s_ms)*np.tan(theta_ms/180.0*np.pi/2.0)*1.1, t_diabolo*2)),
    dl=[2, 2, t_diabolo],
    name="MS"
    # enforce=True
)

refine_ms_wingedge1 = td.MeshOverrideStructure(
    geometry=td.Box(center=(l_ms/2.0, 0.0, 0.0), size=(4.0, (l_ms-s_ms)*np.tan(theta_ms/180.0*np.pi/2.0)*1.1, t_diabolo*5)),
    dl=[0.5, 1.0, t_diabolo],
    name="MS Wing 1"
    # enforce=True
)

refine_ms_wingedge2 = td.MeshOverrideStructure(
    geometry=td.Box(center=(-l_ms/2.0, 0.0, 0.0), size=(4.0, (l_ms-s_ms)*np.tan(theta_ms/180.0*np.pi/2.0)*1.1, t_diabolo*2)),
    dl=[0.5, 1.0, t_diabolo],
    name="MS Wing 2"
    # enforce=True
)


refine_msboxtop = td.MeshOverrideStructure(
    geometry=td.Box(center=(0.0, 0.0, t_magbox), size=(l_magbox*1.1, d_ms*1.1, t_diabolo*2)),
    dl=[1.0, 1.0, t_diabolo],
    name="MS Box Top"
    # enforce=True
)

refine_roi_box = td.MeshOverrideStructure(
    # geometry=td.Box(center=(0.0, 0.0, t_magbox/2.0), size=(l_ms*1.1, (l_ms-s_ms)*np.tan(theta_ms/180.0*np.pi/2.0)*1.1, t_roi)),
    geometry=td.Box(center=(0.0, 0.0, t_magbox/2.0), size=(l_magbox*1.1, d_ms*1.1, t_roi)),
    dl=[0.5, 0.5, 0.5],
    name="ROI_BOX"
    # enforce=True
)

refine_roi_ms = td.MeshOverrideStructure(
    geometry=td.Box(center=(0.0, 0.0, t_magbox/2.0), size=(l_ms*1.1, (l_ms-s_ms)*np.tan(theta_ms/180.0*np.pi/2.0)*1.1, t_roi)),
    # geometry=td.Box(center=(0.0, 0.0, t_magbox/2.0), size=(l_magbox*1.1, d_ms*1.1, t_roi)),
    dl=[5.0, 5.0, 5.0],
    name="ROI_MS"
    # enforce=True
)

refine_rfstrip_edge = td.MeshOverrideStructure(
    geometry=(td.Box(center=(w_rfstrip/2.0, 0.0, t_ufch+t_rfstrip/2.0), size=(t_rfstrip*2, l_ms*0.8, t_rfstrip*2))+
              td.Box(center=(-w_rfstrip/2.0, 0.0, t_ufch+t_rfstrip/2.0), size=(t_rfstrip*2, l_ms*0.8, t_rfstrip*2))),
    dl=[t_rfstrip, t_rfstrip, t_rfstrip],
    # enforce=True
)

refine_rfstrip_edge1 = td.MeshOverrideStructure(
    geometry=td.Box(center=(w_rfstrip/2.0, 0.0, t_ufch+t_rfstrip/2.0), size=(t_rfstrip*4, td.inf, t_rfstrip*4)),
    dl=[t_rfstrip, t_rfstrip, t_rfstrip],
    name="RF Stripline Edge1"
    # enforce=True
)
refine_rfstrip_edge2 = td.MeshOverrideStructure(
    geometry=td.Box(center=(-w_rfstrip/2.0, 0.0, t_ufch+t_rfstrip/2.0), size=(t_rfstrip*4, td.inf, t_rfstrip*4)),
    dl=[t_rfstrip, t_rfstrip, t_rfstrip],
    name="RF Stripline Edge2"
    # enforce=True
)


# define a GridSpec
grid_spec = td.GridSpec.auto(
    min_steps_per_wvl=10, wavelength=lda0, override_structures=[refine_roi_ms, 
                                                                refine_roi_box, 
                                                                refine_ms,  
                                                                refine_ms_wingedge1,
                                                                refine_ms_wingedge2,
                                                                refine_rfstrip_edge1, 
                                                                refine_rfstrip_edge2]
)

run_time = 200.0/fwidth # simulation run time

# define simulation
Lz_struct = t_sub + t_uf
size_sim = (p_ms, p_ms, Lz_struct + 2.2 * offset)
sim = td.Simulation(
    size=size_sim,
    center=(0, 0, Lz_struct/2.0-t_sub),
    grid_spec=grid_spec,
    structures=[substrate, ufluidics, magbox, *msatom, rfstrip],
    sources=[plane_wave],
    monitors=[S11_monitor, S21_monitor, 
              field_monitor_xc0, field_monitor_yc0, 
              field_monitor_zc0, field_monitor_zct, 
              field_monitor_ms_boxcenter, field_monitor_roi],
    run_time=run_time,
    boundary_spec=boundary_spec,
    symmetry=(-1, 1, 0),  # symmetry is used to reduce the computational load
)

sim.plot_3d()

Cross-section of Meshes

f, ax = plt.subplots(1, 3, tight_layout=True, figsize=(15, 5))


sim.plot_grid(z=t_uf+offset, ax=ax[0])

# metasurface layer
sim.plot(z=z_fp_compen, ax=ax[1])
sim.plot_grid(z=z_fp_compen, ax=ax[1])
# ax[1].set_xlim(-d_ms, d_ms)
# ax[1].set_ylim(-d_ms, d_ms)

sim.plot(y=0, ax=ax[2])
sim.plot_grid(y=0, ax=ax[2])
ax[2].set_ylim(-(t_sub+t_uf), t_sub+t_uf)
plt.show()

# job = web.Job(simulation=sim, task_name="THz Slotted Diabolo Stripline")
# job.estimate_cost()

# filename = f"data/simulation_data_{job.task_id}.hdf5"
# sim_data = job.run(path=filename)

# web.real_cost(job.task_id)

S11 and S21

# task_id = job.task_id
task_id = "fdve-1f9be6ba-0e46-4a6f-8a2a-c152b7d54102"
filename = f"data/simulation_data_{task_id}.hdf5"
# sim_data = td.web.load(task_id, filename)
sim_data = td.SimulationData.from_file(fname=filename)
S11 = sim_data["S11"].flux
S21 = sim_data["S21"].flux
# Create a Plotly figure
fig = go.Figure()

# Plot S11
fig.add_trace(go.Scatter(
    x=freqs_GHz,
    y=S11,
    mode='lines',
    name="S11",
    line=dict(color='blue')
))

# Plot S21 in dB
fig.add_trace(go.Scatter(
    x=freqs_GHz,
    y=20 * np.log10(S21),
    mode='lines',
    name="S21",
    line=dict(color='red')
))

# Update the layout with labels, titles, etc.
fig.update_layout(
    title=f"S11 and S21 Spectrum",
    xaxis_title="Frequency [GHz]",
    yaxis_title="S-parameters [dB]",
    legend=dict(title="S-parameters "),
    template='plotly_white',  # You can choose different templates like 'plotly', 'ggplot2', etc.
    width=800,
    height=600
)

# Show the plot
fig.show()

Field Cross Section

f, axs = plt.subplots(2, 3, tight_layout=True, figsize=(15, 10))

valchoice = "abs"

# Transverse Electric field
xrange = l_ms/1.5
yrange = l_ms/1.5
sim_data.plot_field(field_monitor_name="field_xc0", field_name="Ex", val=valchoice, ax=axs[0, 0])
axs[0, 0].set_xlim(-xrange, xrange)
axs[0, 0].set_ylim(-yrange, yrange)
sim_data.plot_field(field_monitor_name="field_yc0", field_name="Ex", val=valchoice, ax=axs[0, 1])
axs[0, 1].set_xlim(-xrange, xrange)
axs[0, 1].set_ylim(-yrange, yrange)
sim_data.plot_field(field_monitor_name="field_zc0", field_name="Ex", val=valchoice, ax=axs[0, 2])

# Transverse Magnetic Field
sim_data.plot_field(field_monitor_name="field_xc0", field_name="Hy", val=valchoice, ax=axs[1, 0])
axs[1, 0].set_xlim(-xrange, xrange)
axs[1, 0].set_ylim(-yrange, yrange)
sim_data.plot_field(field_monitor_name="field_yc0", field_name="Hy", val=valchoice, ax=axs[1, 1])
axs[1, 1].set_xlim(-xrange, xrange)
axs[1, 1].set_ylim(-yrange, yrange)
sim_data.plot_field(field_monitor_name="field_zc0", field_name="Hy", val=valchoice, ax=axs[1, 2])

monitor_name = "field_roi"

x_coords = sim_data[monitor_name].Hy["x"]
y_coords = sim_data[monitor_name].Hy["y"]
z_coords = sim_data[monitor_name].Hy["z"]
f_coords = sim_data[monitor_name].Hy["f"]

f_probe = 463*GHz
idx_f = np.argmin(np.ravel(np.abs(f_coords-f_probe)))

# probe position
x_probe = 0.0
y_probe = 0.0
z_probe = t_magbox/2.0

idx_x = np.argmin(np.ravel(np.abs(x_coords-x_probe)))
idx_y = np.argmin(np.ravel(np.abs(y_coords-y_probe)))
idx_z = np.argmin(np.ravel(np.abs(z_coords-z_probe)))


Hfield_y_xcross = np.abs(sim_data[monitor_name].Hy[:,idx_y,idx_z,idx_f])
Hfield_y_ycross = np.abs(sim_data[monitor_name].Hy[idx_x, :,idx_z,idx_f])
Hfield_y_zcross = np.abs(sim_data[monitor_name].Hy[idx_x,idx_y,:,idx_f])


Hfield_x_probe = np.abs(sim_data[monitor_name].Hx[idx_x,idx_y,idx_z,:])
Hfield_y_probe = np.abs(sim_data[monitor_name].Hy[idx_x,idx_y,idx_z,:])
Hfield_z_probe = np.abs(sim_data[monitor_name].Hz[idx_x,idx_y,idx_z,:])


# probe position
x_probe_corner = 5
y_probe_corner = 5
z_probe_corner = 4

idx_x_corner = np.argmin(np.ravel(np.abs(x_coords-x_probe_corner)))
idx_y_corner = np.argmin(np.ravel(np.abs(y_coords-y_probe_corner)))
idx_z_corner = np.argmin(np.ravel(np.abs(z_coords-z_probe_corner)))


Hfield_y_xcross_corner = np.abs(sim_data[monitor_name].Hy[:,idx_y_corner,idx_z_corner,idx_f])
Hfield_y_ycross_corner = np.abs(sim_data[monitor_name].Hy[idx_x_corner, :,idx_z_corner,idx_f])
Hfield_y_zcross_corner = np.abs(sim_data[monitor_name].Hy[idx_x_corner,idx_y_corner,:,idx_f])


Hfield_x_probe_corner = np.abs(sim_data[monitor_name].Hx[idx_x_corner,idx_y_corner,idx_z,:])
Hfield_y_probe_corner = np.abs(sim_data[monitor_name].Hy[idx_x_corner,idx_y_corner,idx_z,:])
Hfield_z_probe_corner = np.abs(sim_data[monitor_name].Hz[idx_x_corner,idx_y_corner,idx_z,:])

width_beam = 2000.0 #um
normalizefactor = np.sqrt(0.04*p_ms/width_beam)*td.MU_0*1E12*1E4

fig = go.Figure(
    data=[
        go.Scatter(x=x_coords, y=Hfield_y_xcross*normalizefactor, name=f"(y, z)=({y_probe}, {z_probe})"),
        go.Scatter(x=y_coords, y=Hfield_y_ycross*normalizefactor, name=f"(z, x)=({z_probe}, {x_probe})"),
        go.Scatter(x=z_coords, y=Hfield_y_zcross*normalizefactor, name=f"(x, y)=({x_probe}, {y_probe})"), 
        go.Scatter(x=x_coords, y=Hfield_y_xcross_corner*normalizefactor, name=f"(y, z)=({y_probe_corner}, {z_probe_corner})"),
        go.Scatter(x=y_coords, y=Hfield_y_ycross_corner*normalizefactor, name=f"(z, x)=({z_probe_corner}, {x_probe_corner})"),
        go.Scatter(x=z_coords, y=Hfield_y_zcross_corner*normalizefactor, name=f"(x, y)=({x_probe_corner}, {y_probe_corner})"), 
        ],
    layout=dict(
        title=f"Cross section at f={round(float(f_coords[idx_f]/GHz))} GHz",
        xaxis_title="x/y/z coord",
        # yaxis_title="Hy [abs]",
        yaxis_title="By [G]",
        legend_title="Cross Section",
        width=900, 
        height=600,
        xaxis_range=[-5,t_magbox],
        # yaxis_range=[0.0037,0.0040]
        yaxis_range=np.array([.97, 1.03])*float(Hfield_y_zcross[idx_z])*normalizefactor
    )
)
fig.show()

fig = go.Figure(
    data=[
        go.Scatter(x=freqs_GHz, y=Hfield_x_probe*normalizefactor, name=f"Hx"),
        go.Scatter(x=freqs_GHz, y=Hfield_y_probe*normalizefactor, name=f"Hy"),
        go.Scatter(x=freqs_GHz, y=Hfield_z_probe*normalizefactor, name=f"Hz"), 
        ],
    layout=dict(
        title=f"Field Probe at ({x_probe}, {y_probe}, {z_probe})",
        xaxis_title="Frequency [GHz]",
        yaxis_title="By [G]",
        legend_title="Component",
        width=900, 
        height=600,
    )
)
fig.show()

Magnetic Field inside Slot

phase_x = np.angle(sim_data[monitor_name].Hx[:,:,:,idx_f])
phase_y = np.angle(sim_data[monitor_name].Hy[:,:,:,idx_f])
phase_z = np.angle(sim_data[monitor_name].Hz[:,:,:,idx_f])
amp_x = sim_data[monitor_name].Hx[:,:,:,idx_f]
amp_y = sim_data[monitor_name].Hy[:,:,:,idx_f]
amp_z = sim_data[monitor_name].Hz[:,:,:,idx_f]

# probe position
x1 = -4
y1 = -4
z1 = 0
x2 = 4
y2 = 4
z2 = 7.5

idx_x1 = np.argmin(np.ravel(np.abs(x_coords-x1)))
idx_y1 = np.argmin(np.ravel(np.abs(y_coords-y1)))
idx_z1 = np.argmin(np.ravel(np.abs(z_coords-z1)))

idx_x2 = np.argmin(np.ravel(np.abs(x_coords-x2)))+1
idx_y2 = np.argmin(np.ravel(np.abs(y_coords-y2)))+1
idx_z2 = np.argmin(np.ravel(np.abs(z_coords-z2)))+1

x_mesh, y_mesh, z_mesh = np.meshgrid(x_coords[idx_x1:idx_x2], y_coords[idx_y1:idx_y2], z_coords[idx_z1:idx_z2])
x_meshflat = np.ravel(x_mesh)
y_meshflat = np.ravel(y_mesh)
z_meshflat = np.ravel(z_mesh)
shift = np.exp(-1j*np.angle(sim_data[monitor_name].Hy[idx_x,idx_y,idx_z,idx_f]))
Hx_flat = np.ravel(amp_x[idx_x1:idx_x2, idx_y1:idx_y2, idx_z1:idx_z2])*shift
Hy_flat = np.ravel(amp_y[idx_x1:idx_x2, idx_y1:idx_y2, idx_z1:idx_z2])*shift
Hz_flat = np.ravel(amp_z[idx_x1:idx_x2, idx_y1:idx_y2, idx_z1:idx_z2])*shift
# num_frame = 2
# phaseshift = np.linspace(0, 2*np.pi, num_frame, endpoint=False)
# shiftor = np.exp(-phaseshift*1j)
# Hfield_timelapse = [ np.array([Hx_flat, Hy_flat, Hz_flat])*ss for ss in shiftor]

width_beam = 2000.0 #um
normalizefactor = np.sqrt(0.04*p_ms/width_beam)*td.MU_0*1E12*1E4
cenvalue = float(Hfield_y_zcross[idx_z])*normalizefactor
skip = 1
xx=x_meshflat[1:-1:skip]
yy=y_meshflat[1:-1:skip]
zz=z_meshflat[1:-1:skip]
field_x = Hx_flat.real[1:-1:skip]*normalizefactor
field_y = Hy_flat.real[1:-1:skip]*normalizefactor
field_z = Hz_flat.real[1:-1:skip]*normalizefactor

crange = np.array([.98, 1.02])*cenvalue
conesize = float(np.mean(np.sqrt(field_x**2+field_y**2+field_z**2))*len(xx)**(1/3.0)/5.0)
fig = go.Figure(data=go.Cone(
    x=xx,
    y=yy,
    z=zz,
    u=field_x,
    v=field_y,
    w=field_z,
    autocolorscale=True,
    colorscale="Plasma",
    opacity=0.5,
    sizemode="absolute",
    sizeref=conesize,
    colorbar=dict(
        title="B [G]",  # Title for the colorbar
        titleside="top",
        ticks="outside",
    )
))

fig.update_layout(
    width=900, 
    height=600,
    scene=dict(
        aspectratio=dict(x=1, y=1, z=0.8),
        camera_eye=dict(x=1.2, y=1.2, z=0.6),
        xaxis=dict(
            title='X [um]',  # Title for the X-axis
            showbackground=False,
            showgrid=False
        ),
        yaxis=dict(
            title='Y [um]',  # Title for the Y-axis
            showbackground=False,
            showgrid=False
        ),
        zaxis=dict(
            title='Z [um]',  # Title for the Z-axis
            showbackground=False,
            showgrid=False
        )
    )
)

fig.show()

Histogram of Field Distribution

# create the bins
counts, bins = np.histogram(field_y, bins=np.linspace(np.min(field_y), np.max(field_y), 50))
bins = 0.5 * (bins[:-1] + bins[1:])

fig = px.bar(x=bins, y=counts, labels={'x':'B_y [G]', 'y':'count'})

fig.update_layout(
    width=900, 
    height=600)
fig.show()

percent = 1.0
pmpercent_lower = np.mean(field_y)*(1-percent/100.0)
pmpercent_upper = np.mean(field_y)*(1+percent/100.0)
idx_ll = np.argmin(np.abs(bins-pmpercent_lower))
idx_uu = np.argmin(np.abs(bins-pmpercent_upper))
portion_inpm = np.sum(counts[idx_ll:idx_uu])/np.sum(counts)

# fig = go.Figure(
#     data=[go.Cone(
#                 x=x_meshflat[1:-1:4],
#                 y=y_meshflat[1:-1:4],
#                 z=z_meshflat[1:-1:4],
#                 u=Hfield_timelapse[0][0].real[1:-1:4],
#                 v=Hfield_timelapse[0][1].real[1:-1:4],
#                 w=Hfield_timelapse[0][2].real[1:-1:4],
#                 opacity=0.5,
#                 # colorscale='Blues',
#                 colorscale="Portland",
#                 sizemode="absolute",
#                 hoverinfo="text+u+v+w",
#                 text="H Field",
#                 sizeref=float(conesize))],
#     layout=go.Layout(
#         width=900, 
#         height=600,
#         scene=dict(aspectratio=dict(x=1, y=1, z=0.8),
#                              camera_eye=dict(x=1.2, y=1.2, z=0.6)),
#         title_text="H Field Timelapse", hovermode="closest",
#         updatemenus=[dict(type="buttons",
#                           buttons=[dict(label="Play",
#                                         method="animate",
#                                         args=[None])])]
#                                         ),
#     frames=[
#         go.Frame(
#         data=[
#             go.Cone(
#                 x=x_meshflat,
#                 y=y_meshflat,
#                 z=z_meshflat,
#                 u=Hfield_timelapse[k][0].real,
#                 v=Hfield_timelapse[k][1].real,
#                 w=Hfield_timelapse[k][2].real,
#                 opacity=0.5,
#                 # colorscale='Blues',
#                 colorscale="Portland",
#                 sizemode="absolute",
#                 hoverinfo="text+u+v+w",
#                 text="H Field",
#                 sizeref=float(conesize))
#             ])
#         for k in range(num_frame)]
# )


# fig.show()